;+
; NAME:
;       PAPER_PLOT_LV
;
; PURPOSE:
;       Make the L-V plot for the EMAF paper.  Uses 12CO (Dame 2001),
;       13CO (GRS), and HI (IGPS).
;
; CATEGORY:
;       EMAF plotting
;
; CALLING SEQUENCE:
;       PAPER_PLOT_LV
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE  [Creates 3 separate .eps file]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  10/12/11, TPEB -- Initial version.
;       Modified: 11/14/12, TPEB -- Moved to emaf/ directory and made
;                                   compliant.
;       Modified: 01/29/13, TPEB -- Moved more commands to CG library
;                                   routines.
;       Modified: 03/18/13, TPEB -- Added B&W version for published
;                                   paper to avoid paying $$$.
;
;-

PRO PAPER_PLOT_LV
  
  
  ;; Do twice -- first in color, next in BW
  FOR colortable = 41, 42 DO BEGIN
  
  
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  
  IF n_elements(emaf) EQ 0 THEN restore,'./emaf_paper/bgps_emaf.sav',/ver
  
  ;; Parallax Data
  restore,'./ancillary/bessel_parallax.sav',/ver
  
  ;; Info for the 2 loops!
  smol = ['co','hi','13co']
  lvfn = ['./local/co_surveys/COGAL_lv_deep_mom_smoo.fits',$
          './local/hi_lv/HIGAL_LV_plot.fits',$
          './local/13co_lv/13COGAL_LV_plot.fits']
  scaling = ['log   ','','log   ']
  
  IF colortable EQ 41 THEN BEGIN
     fncol = ''
     bw = 0b
     mess = ' in color with CT 41'+string(colortable,format="(I0)")
     masercol = 'Bisque'
     reidcol = 'grn7'
  ENDIF ELSE BEGIN
     fncol = '_bw'
     bw = 0b
     mess = ' in black/white with CT '+string(colortable,format="(I0)")
     masercol = 'blk2'
     reidcol='blk7'
     
  ENDELSE
  
  ;; Loop through surveys
  FOR mol=0,n_elements(smol)-1 DO BEGIN
     
     message,'Now making plot for '+smol[mol]+mess,/inf
     
     
     myps,'emaf_paper/plots/plot_lv_'+smol[mol]+fncol+'.eps',$
          ct=0,xsize=6.5,ysize=3.83,BW=bw
     
     lv = readfits(lvfn[mol],hdr)
     
     IF mol NE 1 THEN lv = alog10(lv)
     CASE mol OF
        0: range = [-1.6,1.1]
        1: range = set_plot_range(lv)
        2: range = [-1.6,0]
     ENDCASE
     vrange = [-20,265]         ;-60
     
     ;; Charsize setting
     cs = 0.7
     
     ;;=====================================================
     ;; Plot the 12CO map
     plotmap,lv,hdr,range=range,/notsquare,axsel=7,ct=colortable,xr=[70,5],$
             yr=vrange,xmargin=[8,12],charsize=cs,min_dpi=300,ymargin=[4,1]
     

     cgLoadct,colortable,/silent
     
     histz = vrange[1] - 70.

     ;; cgColorFill,[20,7,7,20],[histz,histz,histz+50,histz+50],color=1
     ;; cgColorFill,[7.5,6,6,7.5],[histz+15,histz+15,histz+50,histz+50],color=1
     ;; cgColorFill,[35,40,40,35],[135,135,160,160],color=1
     
     ;;=====================================================
     ;; PLOT Rotation Curves
     npts = 1001
     ells = findgen(npts)/10.
     vtan = fltarr(npts)
     dtan = fltarr(npts)
     
     ;; R09 With VS = counter-rotating
     dtan = !MW.R0*cos(ells*!dtor)
     vtan = vphys2vlsr(ells,0.,dtan)
     cgOplot,ells,vtan,thick=5,color=reidcol,linestyle=3  
     print,m4_stat(ells),m4_stat(vtan)
     
     ;; ;; R09 With VS = 0
     ;; dtan = !MW.R0*cos(ells*!dtor)
     ;; vtan = vphys2vlsr(ells,0.,dtan, VS=0.)
     ;; cgOplot,ells,vtan,thick=5,color='GRN7',linestyle=3  
     ;; print,m4_stat(ells),m4_stat(vtan)
     
     ;; Gits & Shiggles -- Clemens
     dtan = !MW.R0*cos(ells*!dtor)
     vtan = vphys2vlsr(ells,0.,dtan, VS=0., /CLEM)
     cgOplot,ells,vtan,thick=5,color='WT1',linestyle=3  
     print,m4_stat(ells),m4_stat(vtan)
     
     
     ;; Indicate W43 Region
     hthick = 3
     cgPlots,[29.5,29.5,31.5,31.5,29.5],[110,80,80,110,110],thick=hthick+2,$
             color='wt1'
     
     
     ;; ;;Plot Near 3-kpc arm
     ;; col3 = 'tg1'
     ;; cgPlots,[5,23],[-32.3,42.58],thick=hthick+2,color=col3,linestyle=0
     ;; cgPlots,[5,23],[-32.3,42.58]+19.7/2,thick=hthick+2,color=col3,linestyle=0
     ;; cgPlots,[5,23],[-32.3,42.58]-19.7/2,thick=hthick+2,color=col3,linestyle=0
     
     
     
     ;;=====================================================
     ;; Hash out the regions associated with the Bar and innards of
     ;;   the Galaxy.
     ;; Upper region
     cgPlots,21,150
     cgPlots,21,85,/cont,thick=hthick,linestyle=0,color='Black'
     cgPlots,7.5,40,/cont,thick=hthick,linestyle=0,color='Black'
     cgPlots,5,31.66666666,/cont,thick=hthick,linestyle=1,color='Black'
     ;; Lower region
     cgPlots,5,-5.55555555
     cgPlots,7.5,0,/cont,thick=hthick,linestyle=1,color='Black'
     cgPlots,21,30,/cont,thick=hthick,linestyle=0,color='Black'
     cgPlots,21,vrange[0],/cont,thick=hthick,linestyle=0,color='Black'
     
     ;; Fill with hashing
     cgColorFill,[21,21,7.5,7.5],[150,85,40,histz],/line_fill,$
              color='Black',orientation=45,thick=hthick

     cgColorFill,[21,21,7.5,7.5],[vrange[0],30,0,vrange[0]],/line_fill,$
              color='Black',orientation=45,thick=hthick
     
     cgLoadct,colortable,/silent
     integral = '!9' + string("362B) + '!X'
     ;;"
     
     cgColorbar,/vertical,charsize=cs,/right,range=range,$
                position=[0.90,!y.window[0],0.93,!y.window[1]],minor=5,$
                title=scaling[mol]+integral+' T!dmb!n db  [K arcdeg]',$
                xthick=3,ythick=3,divis=0
     
     
     
     
     ;; ;;=====================================================
     ;; ;; HI Tangencies
     ;; readcol,'./rot_curve/lbvtan.txt',hil,hib,hiv,format='F,F,F'
     ;; cgPlotsym,8,0.2,/fill
     ;; oplot,hil,hiv,psym=8,color=cgColor('ORG8')
     
     
     
     ;; ;; Attempted Following of contours
     ;; ells = [33,29,27,25,22,19,17,15,13,12,10,9,7,5,2]
     ;; vtan = [107,111,113,116,120,129,139,146,156,164,175,183,191,182,160]
     ;; oplot,ells,vtan,thick=5,color=cgColor('TG2'),linestyle=3
     ;; oplot,ells,vtan+6,thick=5,color=cgColor('GRN7'),linestyle=3
     
     
     ;;=====================================================
     ;; Plot BGPS objects
     cgOplot,emaf.l,emaf.vlsr,psym=16,color='Opposite',symsize=0.3
     
     
     ;;=====================================================
     ;; Plot BeSSeL masers
     name=['G23.0-0.4','G23.4-0.2','G23.6-0.1','W51 IRS2','W51 Main']
     masi = WHERE_ARRAY(name,parallax.name,nmas)
     print,nmas
     cgOplot,parallax[masi].l,parallax[masi].vlsr,psym=46,color=masercol
     cgOplot,parallax[masi].l,parallax[masi].vlsr,psym=45,color='Black'
     
     
     ;;=====================================================
     ;; Plotting the GLON histogram
     lh = histogram(emaf.l, /L64, loc=loc,min=7.0,binsize=1.0)
     n = n_elements(loc)
     lh = [0,lh,0]
     dloc = loc[1]-loc[0]
     loc = [loc[0]-dloc,loc,loc[n-1]+dloc]
     
     print,'HISTOGRAM: ',m4_stat(lh)
     
     
     thick=3
     
     ;; histz = vrange[0] + 10.
     histx = 64.
     cgOplot,loc+0.5,lh/1.+histz,psym=10,thick=thick,color='Opposite'
     cgPlots,histx,histz
     cgPlots,histx-43-14.,histz,/cont,linestyle=0,color='Opposite',$
           thick=thick
     cgPlots,histx,histz
     cgPlots,histx,histz+60,/cont,linestyle=0,color=cgColor('Opposite'),$
           thick=thick
     
     cgText,histx+2.5,histz+30,'N per 1'+cgSymbol('deg')+' bin',$
            ALIGNMENT=0.5,ORIENTATION=90,charsize=cs
     
     
     ;; Tick marks
     FOR j=0L,6 DO BEGIN
        IF (j MOD 2) EQ 0 THEN BEGIN
           cgText,histx+0.4,$
                  histz+10*j-3,string(10*j,format="(I0)"),charsize=0.9*cs,$
                  align=1,color=cgColor('Opposite')
           tlen = 1.0
        ENDIF ELSE tlen = 0.5
        cgPlots,histx,histz+10*j
        cgPlots,histx-tlen,histz+10*j,/cont,linestyle=0,color='Opposite',thick=3
        
     ENDFOR
     
     ;; ;; Plot outer spiral arms
     ;; cgPlots,19,0
     ;; cgPlots,70,-70,/cont,linestyle=2
     ;; cgPlots,0,0
     ;; cgPlots,60,-95,/cont,linestyle=2
     
     myps,/done
     
  ENDFOR
  
  cgLoadct,0,/silent
  
  ENDFOR
  
END
